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Turbulence modeling for separated flow 

By P. A, Durbin 

1. Motivation and objectives 

Two projects are described in this report. The first involves assessing turbulence 
models in separated flow. The second addresses the anomalous behavior of certain 
turbulence models, such as k - e, in stagnation point flow. 

The primary motivation for developing turbulent transport models is to provide 
tools for computing non-equilibrium, or complex, turbulent flows. Simple flows 
can be analyzed using data correlations or algebraic eddy viscosities, but in more 
complicated flows such as a massively separated boundary layer, a more elaborate 
level of modeling is required. It is widely believed that at least a two-equation 
transport model is required in such cases. The transport equations determine the 
evolution of suitable velocity and time-scales of the turbulence. The appropriate 
velocity scale for turbulent transport toward a wall is the normal component, v 2 , 
not the turbulent intensity, k. This and other considerations motivated the k — 
e — v 2 model, which can be used in wall-bounded flows. Applications to complex 
geometry require a generalized interpretation of the velocity scale v 2 . The model 
originally was developed for attached or mildly separated boundary layers. Rather 
promising results were obtained in tests of the formulation. Here the model is 
assessed in strongly separated flows. Doing so requires that v 2 be regarded simply as 
a velocity scale that satisfies boundary conditions suitable for the normal component 
of turbulent intensity; it cannot be regarded as the ‘y-component’ because that 
would be inappropriate in general geometries that can have surfaces aligned in any 
direction. This loosened understanding of v 2 presents no operational difficulties. 

Two equation models as well as the k - e - v 2 model predict an anomalously 
large growth of turbulent kinetic energy in stagnation point flows. Even when the 
stagnation point region is not of interest per se , this spurious behavior can upset 
the rest of the flow computation. A formal upper limit to the turbulent time-scale 
in such models alleviates their stagnation point anomaly. This bound is derived 
in §2.4 and is illustrated with the stagnation point flow at the leading edge of an 
airfoil. 

2. Accomplishments 

The governing equations of the k — e — v 2 model will not be presented here. Their 
initial development is described in Durbin (1991), and a more complete description 
of the present work will appear in Durbin (1994). The mean flow satisfies the 
incompressible Navier-Stokes equations with _an eddy viscosity. The turbulence 
model uses the standard k — £ equations, a v 2 transport equation and an elliptic 
relaxation equation for the source term fii in the v 2 -equation. 

The computations were done with an extended version of the INS2D code of 
Rogers and Kwak (1990). The extensions required for the present computations 
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FIGURE 1. Skin friction coefficient on wall downstream of the backstep compared 

to experiments of Jovic and Driver ( , * ) and Driver and Seegmiller ( , 

• ). line is solution to SSG model for Jovic and Driver flow. 

consist of providing subroutines to solve the transport and elliptic relaxation equa- 
tions of the A: — e — v 2 model. This code development will not be described here. 
The numerics are discussed in Durbin (1994). The program can solve full second 
moment closure models as well as eddy viscosity models. 

The present study included assessment of second-moment closures in separated 
flow. The results were generally less satisfactory than Ar — e — v 2 . These second- 
moment computations will not be presented. 

Sections 2. 1-2.3 present three separated flows: these illustrate sharp edge separa- 
tion; smooth wall, pressure driven separation; and unsteady vortex shedding. This 
last case was studied at the suggestion of colleagues at Pratt & Whitney. 

2.1 Flow over backward facing steps 

The backstep flow configurations studied experimentally by Jovic and Driver 
(1993) and by Driver and Seegmiller (1985) were computed with the k — e — v 2 model 
and with second-moment closure models. The JD experiment had a step height 
Reynolds number of 5,000; the DS case had Re = 37,500. Both the IP and SSG 
second-moment closures as incorporated into the elliptic relaxation procedure were 
tried. Computations showed that they significantly under-predict the magnitude of 
the reversed flow downstream of the step, as suggested by the dotted skin friction 
curve in Fig. 1: this should be compared to the solid line and triangles. 

Computed and experimental skin friction coefficients on the wall downstream of 
the step are compared in Fig. 1. The computed reattachment point at x = 6.2 step 
heights is in agreement with the data. The relatively large negative skin friction in 
the JD experiment is due to low Reynolds number. The k - e — v 2 model correctly 
shows this sensitivity to Reynolds number. Overall, the agreement with experiment 
is better than has been found using the standard k — e model with wall functions 
(Driver and Seegmiller, 1985); numerous independent computations have shown 
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FIGURE 2. Mean velocity profiles for the Jovic and Driver experiment, 
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FIGURE 3. Mean velocity profiles for the Driver and Seegmiller experiment. 
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that model under-predicts the reattachment length. 

Figs. 2 and 3 show profiles of the U-component of velocity. These profiles are 
plotted in the form 10Z7 + x to display their evolution down the duct. The agreement 
between model and experiment in the JD case is quite good. In the higher Reynolds 
number DS case, the agreement is good for x < 8, but the model solution for the 
boundary layer downstream of reattachment recovers more slowly than the data. 
This slow recovery downstream of reattachment is a universal problem of turbulence 
models shown by Reynolds stress as well as k — e models. 

2,2 A separated diffuser 

Obi et al. (1993) measured the flow in a one-sided, 10° plane diffuser. The 
expansion ratio of 4.7 was sufficient to produce a separation bubble on the sloping 
wall; hence, this provides a test case for smooth, adverse pressure driven separation. 
The entrance to the diffuser consisted of a long plane channel (aspect ratio 35) in 
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FIGURE 4. Mean velocity profiles in the Obi, et al. diffuser. Both Eq. (9) ( ) 

and a constant value of 1.55 ( ) were used for C\ x . The light dashed lines 

show the diffuser surface. 

order to produce fully developed flow. The Reynolds number based on the half- 
height of this channel was 10,000. The computational inflow profiles were obtained 
by solving fully developed channel flow with a parabolic code. Fig. 4 shows profiles 
of 10U + x. The boundary of the diffuser is indicated in the figure— note, however, 
that the aspect ratio of this figure is not unity: the actual duct is more elongated. 
The profiles on the ramp are in good agreement with the data, showing the smooth 
separation; further downstream, the predicted backflow is less strong than the data. 
Note that the mass flux is constant, so less backflow near the lower wall is necessarily 
accompanied by less forward flow in the upper part of the channel. Second moment 
closure computations of this flow failed to predict the separation in this flow. Both 
SSG and IP second-moment closure models gave only a tiny region of reversed flow 
at the foot of the ramp. 

2.S Vortex shedding behind a triangular cylinder 

The flow around a triangular cylinder in a duct was measured by Sjunnesson 
(presented in Johansson et al. 1993). This geometry provides an example of bluff 
body flow with fixed separation points. The study by Sjunnesson was motivated 
by the application to flame holders. The geometry consists of a 6:1 aspect ratio 
equilateral triangular cylinder centered symmetrically in a duct three cylinder sides 
high. The Reynolds number based on the cylinder side and inlet velocity was about 
42,500. 

Both steady state and (statistically) unsteady solutions were computed. By tak- 
ing a large artificial time-step on the order of the shedding period in length, we ob- 
tained symmetric, steady solutions. Upon introducing an asymmetric disturbance 
and taking a smaller step, an oscillatory solution was obtained. Fig. 5 is a com- 
posite showing the time-averaged U-contours of the unsteady computation in the 
upper half and the steady state solution in the lower. It shows that the steady-state 




Separated flow computations 


101 



FIGURE 5. U-contours for a steady calculation (lower half) and time-average of 
an unsteady computation with vortex shedding (upper half). 



Figure 6. Inst ant aneous vorticity contours showing shedding in the time accurate 
computation. 


separation bubble is over twice as long as the averaged, unsteady bubble. 

Instantaneous vorticity contours, in Fig. 6, show the asymmetric flow in the time 
accurate computation. The interaction between the vortex street and the boundary 
layers in the side walls results in eruption of secondary vorticity. The complexity of 
such interactions underscores the need to resolve the coherent vortices in this type 
of flow. 

Fig. 7 shows profiles of the time-averaged U-velocity component in the wake. The 
computational results were obtained by averaging the numerical solution over one 
period of the vortex shedding. The profiles are displayed evenly spaced, but the 
actual locations were x = 0.375, 0.95, 1.525, 3.75, and 9.4 heights downstream 
of the rear face of the cylinder. The agreement between experiment and model is 
excellent. The first profile at x = 0.375 shows the sharp boundary and large velocity 
deficit of the near wake. By x = 0.95 the wake profile has altered substantially 
and is undergoing transition to a Gaussian form. The different curves show the 
convergence of the solution with grid refinement. It can be concluded that the 
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FIGURE 7. Time-averaged velocity profiles in the wake of a triangular cylinder. 

The calculations are shown for three grids: 141 x 101 ( ) 121 x 91 ( )• 

71 x 51 ( ) ’ 


statistical unsteadiness produced by vortex shedding must be resolved in order 
to compute this flow. The steady state computation of this flow is as dubious 
as a homogeneous computation of a spatially homogeneous flow; it leads to quite 
erroneous predictions. 


2.4 On the k — e stagnation point anomaly 

Two equation models as well as k - e - v 2 predict an anomalously large growth 
of turbulent kinetic energy near to stagnation points (Launder and Kato, 1993; 
Menter, 1992). This can cause difficulties in aerodynamic flows with a free-stream 
impinging on a blunt leading edge. Even when the stagnation point region is not of 
interest per ae, this spurious behavior can upset the rest of the flow computation. 
The usual explanation for the stagnation point anomaly is that the eddy viscosity 
formula 


UiUj — 2v t Sij T -k6ij 
o 


( 1 ) 


gives an erroneous normal stress difference (Launder and Kato 1993). In (1), S tj = 
(diUj + djUi )/ 2 is the rate of strain and 


u t = CpkT 


( 2 ) 


is the eddy viscosity. T is the turbulent time-scale (e.g., Jfc/e). 

Some of our computations suggest an alternative understanding of the anomaly: 
as the stagnation point is approached, T becomes very large. The e-equation is of 
the form 


d t e + U • Ve = 



T 


( 3 ) 
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where the rate of turbulent energy production is V — 2vtS x jSji. A large value of 
T in Eq. (3) causes the production of e to be too small, allowing spuriously high 
turbulent kinetic energy. The stagnation point anomaly can be ameliorated by 
imposing a bound on the time-scale. In the following we will derive the constraint 



where SijSji = |S| 2 . In most situations this reduces to T = fc/e; in highly strained 

flow the upper bound comes into play. 

The constraint will be derived by requiring that (1) satisfy 2 k > u 2 > 0, which 
can be called a ‘realizability condition’. The rate of strain tensor Sij is symmetric 
and becomes purely diagonal in principal- axes coordinates. The diagonal elements, 
A q , a = 1 ... 3, are its eigenvalues and satisfy 


A? + A 2 + A 2 = |S| 2 . 

(5) 

In incompressible flow 

( 6 ) 

Ai + A 2 + A 3 = 0 

It follows from (5) and ( 6 ) that 


|A q | = a/|S| 2 /2 

( 7 ) 

in two dimensions (i.e., when A 3 = 0 ), and 


|A«| < \/2|S| 2 /3 

( 8 ) 


in three dimensions. 

If ( 1 ) is written in the principal axes of Sij, it becomes 


— — 2vt\ a + ■ 


h 

3 


( 9 ) 


Of the constrains u 2 a > 0 and 2k > u 2 , Vo, the former is more 
constraint is 2 

2u t max A a < -k . 

O 6 

Substituting (2) into (11) results in the time-scale bound 


stringent; this 


( 10 ) 


which gives 



1 

max X Q 


T < 


2 1 


(ID 


3C M V /W 


( 12 ) 
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Figure 8. Contours of constant (a), with (6) imposed; (b), without 

constraint. Contour intervals of 1.5 x 10~ 3 . 

in two dimensions, and 

T < 

in three dimensions. These bounds might be imposed computationally by Eq. 4. 

Fig. 8 shows k contours for the flow over a NACA4412 airfoil at zero angle of 
attack and with k — 4 x 10 -4 in the approach flow, with and without the con- 
straint on T . This computation was done with the k-e-v 2 model. The constraint 
prevents the spurious growth of k although some amplification still occurs, 

3, Future plans 

Flows with mean swirl are of interest for their role in enhancing mixing both by 
turbulent and mean motion. The swirl can have a stabilizing as well as a destabi- 
lizing effect on the turbulence. I have written an axi-symmetric extension to the 
INS-2D computer program and added the capability to compute swirling flow. This 
is in order to study confined coaxial jets with or without swirl. High swirl can 
produce vortex breakdown on the centerline of the jet. This type of flow occurs in 
various combustors. 

The present vortex shedding calculations suggest that the application of turbu- 
lence models to separation control by external periodic excitation should be ex- 
plored. This is a problem that Hans Kaltenbach has been investigating by LES. 
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